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Abstract - We reconsider the long-standing question of the critical defect hopping rate r^ in the 
one-dimensional totally asymmetric exclusion process (TASEP) with a slow bond (defect). For 
r < rc a. phase separated state is observed due to queuing at the defect site whereas for r > rc the 
defect site has only local effects on the stationary state of the homogeneous system. Mean-field 
theory predicts rc = 1 (when hopping rates outside the defect bond are equal to 1) but numerical 
investigations seem to indicate r^ ~ 0.80(2). Here we improve the numerics to show that r^ > 0.99 
and give strong evidence that indeed Cc = 1 as predicted by mean-field theory, and anticipated by 
recent theoretical findings. 


Introduction. — Despite much progress in recent 
years, our understanding of nonequilibrium stationary 
states is far from complete. This especially concerns the 
effects of disorder and defects in driven diffusive systems. 
Although it is well established that in driven systems al¬ 
ready a localized defect can have a global influence on 
the behavior, many open questions remain. Since e.g. no 
analogue of the Harris criterion [I] is known for nonequi¬ 
librium situations no general statements on the influence 
of weak disorder on the critical behavior can be made [2] . 

Surprisingly even for the simplest model of driven diffu¬ 
sion, the totally asymmetric exclusion process (TASEP), 
the precise influence of a single defect has not been fully 
clarified for a long time. It is well-known since the seminal 
work of Janowsky and Lebowitz mm that such a defect 
has not just a local effect, but changes the nature of the 
stationary state dramatically (for reviews, see e.g. [5j[6]). 
What is not clear up to now is whether a finite critical 
strength of the defect is required to create global effects. 
Mean-held theory predicts a global influence already for 
arbitrarily small defect strengths whereas the most accu¬ 
rate numerical investigations up to date [7] show strong 
indications that a hnite defect strength is required. 

Recently the problem has been newly addressed in the 
mathematical literature, in a series of works [SHU- In 
[8], based on analytical arguments from series expansions 
and results for related systems (e.g. directed polymers) it 


was argued that an arbitrarily small defect in a TASEP 
with open boundaries will have global effects, e.g. on the 
current and on the density prohle. In [9], the authors 
claim to have proved rigorously, that the steady current 
in TASEP with a slow bond is always affected for any 
nonzero defect strength. 

In view of these new hndings, the numerical studies pre¬ 
dicting hnite critical blockage strength, appear even more 
paradoxical, with the = 1 problem hnally settled. It 
remains to understand if the effects of the slow bond are 
so weak that they cannot be observed in numerics, which 
would make the beautiful theoretical result a pure theo¬ 
retical construction without any practical content. 

It is the purpose of this letter to show that this is not 
the case, and to provide detailed results from highly accu¬ 
rate Monte Carlo simulations which strongly support the 
theory. Moreover, we also resolve the apparent numerical 
paradox, revisiting previous numerical studies and point¬ 
ing out exactly where the error of the previous numerical 
studies was. 

Model. — The totally asymmetric simple exclusion 
process (TASEP) is a paradigmatic model of nonequilib¬ 
rium physics (for reviews, see e.g. [T2HT8] ) describes inter¬ 
acting (biased) random walks on a discrete lattice of N 
sites, where an exclusion rule forbids occupation of a site 
by more than one particle. A particle at site k moves to 
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Fig. 1: Open TASEP with a slow bond (r < 1) in the 
middle. 


site fc + 1 at rate p if site fc + 1 is not occupied by another 
particle. The boundary sites k = 1 and k = N are cou¬ 
pled to particle reservoirs. If site 1 is empty, a particle is 
inserted at rate a. A particle on site N is removed from 
the system at rate p. Sites are updated using random- 
sequential dynamics. Throughout the paper, we will set 
p=l. 

Here we consider a system of two TASEPs of length N/2 
coupled by a slow bond between sites N/2 and A^/2-|-l with 
reduced hopping rate r < p (Fig. [1]). This is equivalent to 
a TASEP of N sites with a defect site in the middle and 
defect strength p — r. 

For periodic boundary conditions this problem has been 
analyzed by Janowsky and Lebowitz [314]. Below a crit¬ 
ical rate Tc they found a phase separation into high and 
low density regions due to queuing at the defect site. The 
two phases are separated by a shock (domain wall). The 
phase separation is also reflected in the current-density 
relation (fundamental diagram) which shows a density- 
independent current at intermediate densities due to the 
current-limiting effect of the slow bond. A mean-field the¬ 
ory that neglects correlations at the defect site predicts 
that Tc = 1 13, he. an arbitrarily small defect leads to a 
phase separated stationary state. This is supported by se¬ 
ries expansions obtained from exact results for small sys¬ 
tems 111- Exact results have been obtained for the case 
of sublattice-parallel update with deterministic bulk hop¬ 
ping by Bethe Ansatz m and matrix-product Ansatz [20] . 
Also in this case Tc = 1. 

For open boundary conditions, a mean-field treatment 
of the TASEP with a defect in the middle of the system 
yields Vc = 1 El]. Later Ha et al. [7] studied the problem 
numerically (with rates p = 1, a = P = 1/2). They 
suggested that Vc = 0.80 (2), see next section for details. 

Due to its relevance e.g. for intracellular transport, re¬ 
cently the open TASEP with several defects has attracted 
some attention, see e.g. [22U33 - 

Simulations. — In order to analyze the system with 
a defect we perform Monte Carlo (MC) simulations for 
an open TASEP on a chain of large size {N < 200.000), 
to minimize finite-size effects as much as possible. Both 
bulk hopping rates and boundary hopping rates are cho¬ 
sen equal to 1. This corresponds to a point in the max¬ 
imal current phase which is characterized by a spatially 
homogeneous steady state with bulk density of particles 
Pbuik = 1/2 and the current j(oo) = Pbuik(l - Pbuik) = 1/4 
in the limit of an infinite system size |13| . 



Fig. 2: Finite-size corrections of the steady current for a 
homogeneous model without slow bond. Error bars show 
the 99% conhdence bound, the red line marks the exact 
leading finite-size correction in 1/A^. 


Initially N/2 particles are distributed randomly within a 
homogeneous chain without a defect. Then, the relaxation 
of the system is performed for lOOTV^ single Monte Carlo 
updates [23 , according to the dynamical rules for a = /? = 

1 . 

After the initial relaxation, the weak bond is introduced 
in the middle of the system, meaning that a particle hop 
from site k = N/2 to A^/2 -|- 1 occurs with reduced rate 
r < 1, see Fig. [T] Then the system is relaxed further 
for lOOA^^ single Monte Carlo updates and the average 
current is recorded. It can be written as 

jFs{r,N,a,P) =j{r,oo) + Sps{r,N,a,P) (1) 

where (5 fs are finite-size corrections. The relaxation to 
steady state is controlled by a comparison of the finite- 
size corrections of the current measured numerically with 
the theoretically-predicted value, see Fig. El 

We aim at determining a lower bound for the critical 
Tc at which the phase separation starts. It is well-known 
[35] that the leading finite-size corrections to the current 
of the homogeneous TASEP for a = /3 = 1 are positive, 

jpsir = l,N,a = l,P = 1) = j {r = l,oo) +-^+0{N‘^) 

> j(r = I,oo). (2) 

Therefore, if for some defect hopping rate tq the steady 
current within the error bars is smaller than its limiting 
value, j{ro,N) < 1/4, this would definitely mean that 
fc > cq. The lower bound r* for is then calculated as 
the point where j{r*,N) = 1/4, accounting also for error 
bars, see Eig. [3 Note that our reasoning does not involve 
any a priori assumptions except from the positiveness of 
finite-size corrections to the current. In this way, we obtain 
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Fig. 3: Finite-size current jps (see eq. ([T])) versus defect 
rate r for a = ^ = 1 . With j (r) < j ps i^) the plot 
shows Tc > 0.86 which is significantly larger than Vc = 
0.80(2) suggested in [7]. Currents are averaged over all 
sites with 7.5 • lO"^ — 5 • 10® histories. Error bars show the 
99% confidence interval. 


a lower bound for r,, which depends on the system size. For 
larger system size, finite-size corrections are smaller and 
better estimates for the lower bound can be made, see 
Fig. [31 For system size N = 2 - 10® we obtain the following 
lower bound estimate. 


Vc > 0 . 86 . 


(3) 


It is difficult to improve the lower bound (|3|) by a further 
increase of the system size N, because much larger system 
sizes are not numerically accessible. However, already the 
lower bound value 0.86 dehnitely contradicts the result 
of Ha et al. [7], where Vc = 0.80(2) was found. Their 
estimate was based on a different argument. In order to 
understand the reason for the contradiction, we repeated 
Monte Carlo simulations with the same parameters as in 
[7], and in particular using the much smaller system size 
of = 4100. 

The key quantity analyzed by Ha et ah, is defined as 


Ab oc 


1 

Pbulk,segment 


or 

Afc = 2y/jir = l,N)-jir,N). 
It is assumed to obey the scaling form 


(4) 

(5) 


Ab (rc -r) ^ 


( 6 ) 


near the phase boundary. Then, the best straight line 
fit on the double logarithmic plot of A;, versus Cc — r, 
has lead to the conclusion = 0.80(2). We repeated 
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Fig. 4: Double logarithmic plot of Ab versus (cc — r). Pa¬ 
rameters are N = 4100, a = {3 = 0.5. The box indicates 
the regime plotted in Fig. 4a of [7]. Averaging over all 
sites of the system and over 3 • 10® histories is performed. 
Statistical errors are smaller than the symbol size. 


the relevant Monte Carlo simulations for the parameters 
chosen in [7]. Fig. |4| shows the double logarithmic plot of 
Ab versus — r. The square window corresponds to the 
area shown in the original paper, see Fig. 4a of [7]. We 
can see that what might look as a straight line inside the 
window, certainly fails to straighten outside the window. 
Thus the conclusion of [7] of a phase transition at = 
0.80 (2) is not justified. 

It is instructive at this point to stress the importance 
of a choice of an adequate random generator to perform 
the Monte Carlo update. This choice is crucial for pro¬ 
ducing high quality Monte Carlo data [55]. In Fig. (5] 
Monte Carlo data for the current, produced by different 
random number generators, are compared and show sys¬ 
tematic differences. For our simulations throughout this 
paper we are using Mersenne-Twister-generator which is 
known for producing high quality pseudo-random num¬ 
bers. In Fig. |S] it is seen that using the most common 
Park and Miller new minimal standard linear congruen- 
tial generator m leads to a systematic overestimation of 
the current, and might consequently lead to wrong con¬ 
clusions in the subtle TASEP blockage problem. 

We also note, that the point a = /3 = l/2as chosen 
for studying the blockage problem in |7] lies exactly at 
the phase boundary between the maximal current, low 
density and high density phases [I2UI8] . This may lead to 
further complications and an additional reduction of the 
steady current due to fluctuations. We stress that for our 
study we choose a = /3 = 1, i.e. a point well inside the 
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Fig. 5: Finite-size corrections of the current eq. ([T]), for 
a different choices of a random nnmber generators. Pa¬ 
rameters: N = 4100. Triangles correspond to Mersenne 
Twister random number generator and crosses correspond 
to a standard Park and Miller new minimal standard 
(LC-gen, linear congruential) random number generator 
m- Averaging of the current over all bonds and over 
3 • 10® — 5 • 10® histories is performed. Statistical errors are 
smaller than the symbol size. 

maximal current phase far from the boundaries with the 
low density TASEP phase {a = 1/2) and the high density 
TASEP phase (/3 = 1/2). 

Effects of a defect in finite systems: parallel evo¬ 
lution. — As is already mentioned, both mean-field the¬ 
ory and series expansions arguments hint at = 1. As¬ 
suming the existence of an essential singularity at = 1, 
/(I) — j{r) ^ exp(—a/(1 — r)) [8], further improvement of 
the lower bound for the weak bond problem by an increase 
of the system size is a hopeless enterprise: e.g. a numer¬ 
ical proof Tc > 0.9, rc > 0.95, rc > 0.99) with the direct 
method (see Fig. [S]) would require N > 10^°, N > 10^^, 
N > 10^^^ respectively. 

Instead of increasing the system size, we address the 
problem of a critical blockage strength in different way 
by measuring how a TASEP responds to a slow bond, as 
discussed below. 

After the pre-relaxation performed on the homogeneous 
system as described above, we make two copies of the sys¬ 
tem configuration. Then a slow bond is introduced in 
one copy whereas the other remains homogeneous. Both 
copies evolve in time according to the same protocol, i.e. 
using the same set of random numbers for both systems. 
The averaged density profiles of both copies are compared 
after sufficiently large relaxation time. Due to the de¬ 
fect site, one expects a density gradient forming locally in 
the vicinity of the blockage for any r < 1. If the distur¬ 
bance remains local, the state of the system far from the 


blockage will not change, with respect to a homogeneous 
system. In contrast, a non-local disturbance spreading to 
the whole system would lead to a reduction of the global 
current and to phase separation. Thus the current and 
the density profile are sensitive probes for the effects of 
the slow bond and for the possible occurrence of a phase 
separated state. 

Performing extensive MC simulations we are able to see 
a non-local effect of the blockage up to r = 0.99, both in 
steady current and in local particle density far away from 
the blockage, see Fig. [6l Consequently, a presence of the 
weak bond has a small but systematic effect on both the 
current and the local density Uk far away from the blockage 
site. This shows that the blockage produces perturbations 
which do not remain local, but spread over the bulk. 

Conclusion. — We have revisited a long-standing 
problem of a TASEP with a weak bond. Since the ef¬ 
fects are small for small defect strength, this is a subtle 
problem that requires high numerical accuracy in simu¬ 
lations. We have shown that here even the choice of the 
random number generator is crucial. By Monte Carlo sim¬ 
ulations performed on large systems (up to 200.000 sites), 
we have established a new lower bound for the critical de¬ 
fect strength leading to phase separation. For the current, 
we could clearly show a reduction compared to the value 
j(r- = 1) = 1/4 of the homogeneous system for defect hop¬ 
ping rates up to r = 0.86. Therefore we conclude that 
rc > 0.86. Studying a parallel evolution of two initially 
identical systems, with and without a weak bond, we find 
systematic global effects induced by the weak bond for de¬ 
fect hopping rate up to r < 0.99. Our study supports 
the hypothesis that the critical blockage hopping rate is 
Tc = 1, and definitely rules out a previously obtained crit¬ 
ical value Tc = 0.80(2). This indicates that the mean- 
field theory prediction is indeed correct which is impor¬ 
tant since it is used quite frequently also in more complex 
situations, like flows on networks [38ll42] or as effective 
models for highway traffic near ramps [43] . 

* * * 
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